
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <omp.h>
#include <mpi.h>

#define ceild(n,d)  ceil(((double)(n))/((double)(d)))
#define max(x,y)    ((x) > (y)? (x) : (y))
#define min(x,y)    ((x) < (y)? (x) : (y))
#define myabs(x,y)  (((x) > (y))? ((x)-(y)) : ((y)-(x)))

#if !defined(point) 
#define point 7
#endif

#if point == 7
#define  kernel(A)  A[(t+1)%2][z][y][x] = 0.64 * (A[t%2][z][y][x]) + \
										 0.06 * (A[t%2][z][y][x - 1] + A[t%2][z][y - 1][x] + \
												 A[t%2][z - 1][y][x] + A[t%2][z][y][x + 1] + \
												 A[t%2][z][y + 1][x] + A[t%2][z + 1][y][x])
#define XSLOPE 1
#define YSLOPE 1
#define ZSLOPE 1
#elif point == 19
#define  kernel(A) A[(t+1)%2][z][y][x] = 2.666f * A[(t)%2][z][y][x] - \
		0.166f * (A[(t)%2][z-1][y][x] + A[(t)%2][z+1][y][x] + A[(t)%2][z][y-1][x] + A[(t)%2][z][y+1][x] + A[(t)%2][z][y][x+1] + A[(t)%2][z][y][x-1]) - \
		0.0833f * (A[(t)%2][z-1][y-1][x] + A[(t)%2][z+1][y-1][x] + A[(t)%2][z-1][y+1][x] + A[(t)%2][z+1][y+1][x] + \
			   A[(t)%2][z-1][y][x-1] + A[(t)%2][z+1][y][x-1] + A[(t)%2][z][y-1][x-1] + A[(t)%2][z][y+1][x-1] + \
			   A[(t)%2][z-1][y][x+1] + A[(t)%2][z+1][y][x+1] + A[(t)%2][z][y-1][x+1] + A[(t)%2][z][y+1][x+1])
#define XSLOPE 1
#define YSLOPE 1
#define ZSLOPE 1
#elif point == 27
#define kernel(A) A[(t+1)%2][z][y][x] = 0.54 * (A[(t)%2][z][y][x]) + \
										0.03 * (A[(t)%2][z][y][x-1] + A[(t)%2][z][y-1][x] + \
												A[(t)%2][z][y+1][x] + A[(t)%2][z][y][x+1] + \
												A[(t)%2][z-1][y][x] + A[(t)%2][z+1][y][x])+ \
										0.01 * (A[(t)%2][z-1][y][x-1] + A[(t)%2][z-1][y-1][x] + \
												A[(t)%2][z-1][y+1][x] + A[(t)%2][z-1][y][x+1] + \
												A[(t)%2][z][y-1][x-1] + A[(t)%2][z][y+1][x-1] + \
												A[(t)%2][z][y-1][x+1] + A[(t)%2][z][y+1][x+1] + \
												A[(t)%2][z+1][y][x-1] + A[(t)%2][z+1][y-1][x] + \
												A[(t)%2][z+1][y+1][x] + A[(t)%2][z+1][y][x+1])+ \
										0.02 * (A[(t)%2][z-1][y-1][x-1] + A[(t)%2][z-1][y+1][x-1] + \
												A[(t)%2][z-1][y-1][x+1] + A[(t)%2][z-1][y+1][x+1] + \
												A[(t)%2][z+1][y-1][x-1] + A[(t)%2][z+1][y+1][x-1] + \
												A[(t)%2][z+1][y-1][x+1] + A[(t)%2][z+1][y+1][x+1]);
#define XSLOPE 1
#define YSLOPE 1
#define ZSLOPE 1
#endif
#define DATA_TYPE double

// #define CHECK
#ifdef CHECK
#define TOLERANCE  0
#endif


void init(DATA_TYPE *A, int domain_z, int domain_y, int domain_x, int z_start, 
        int z_end, int y_start, int y_end, int x_start, int x_end) {
    for (int k = z_start; k < z_end; k++) {
        for (int j = y_start; j < y_end; j++) {
            for (int i = x_start; i < x_end; i++) {
                A[k * (domain_x * domain_y) + j * domain_x + i] = (DATA_TYPE) (1.0 * (rand() % 1024));
                A[domain_z * domain_y * domain_x + k * (domain_y * domain_x) + j * domain_x + i] = 0;
            }
        }
    }
}

int pack(DATA_TYPE *A, int domain_z, int domain_y, int domain_x, int z_start, 
        int z_end, int y_start, int y_end, int x_start, int x_end, double *buffer, int t) {
    int count = 0;
    for (int k = z_start; k < z_end; k++) {
        for (int j = y_start; j < y_end; j++) {
            for (int i = x_start; i < x_end; i++) {
                buffer[count] = A[(t % 2) * domain_z * domain_y * domain_x + k * (domain_x * domain_y) + j * domain_x + i];
                count++;
            }
        }
    }
    return count;
}

void unpack(DATA_TYPE *A, int domain_z, int domain_y, int domain_x, int z_start, 
        int z_end, int y_start, int y_end, int x_start, int x_end, double *buffer, int t) {
    int count = 0;
    for (int k = z_start; k < z_end; k++) {
        for (int j = y_start; j < y_end; j++) {
            for (int i = x_start; i < x_end; i++) {
                A[(t % 2) * domain_z * domain_y * domain_x + k * (domain_x * domain_y) + j * domain_x + i] = buffer[count];
                count++;
            }
        }
    }
}

int main(int argc, char * argv[]){

	int iam, np;
	double begin, elaspe;
	
	MPI_Status status;
	MPI_Init(&argc,&argv);
	MPI_Comm_rank(MPI_COMM_WORLD,&iam);	
	MPI_Comm_size(MPI_COMM_WORLD,&np);//np=8

	long int t, i, j, k;
	int NX = atoi(argv[1]);
	int NY = atoi(argv[2]);
	int NZ = atoi(argv[3]);
	int T  = atoi(argv[4]);
	int ghost_x = atoi(argv[5]);
	int ghost_y = atoi(argv[6]);
	int ghost_z = atoi(argv[7]);

	int npx = pow(np, 1.0 / 3);
	int npy = npx;
	int npz = npx;
	int ix = ceild(NX,npx);
	int iy = ceild(NY,npy);
	int iz = ceild(NZ,npz);
	
	int domain_x = ix + 2*ghost_x;
	int domain_y = iy + 2*ghost_y;
	int domain_z = iz + 2*ghost_z;

    int tidx = (iam % (npy * npx)) % npx;
    int tidy = (iam % (npy * npx)) / npx;
    int tidz = iam / (npy * npx);
    
    if (iam < npx * npy * npz) {
        if(tidx == npx - 1){ //右边界进程x大小修正
            domain_x = NX - (ix * (npx-1)) + 2*ghost_x;
        }
        if(tidy == npy - 1){ //上边界进程y大小修正
            domain_y = NY - (iy * (npy-1)) + 2*ghost_y;
        }
        if(tidz == npz - 1){ //顶边界进程z大小修正
            domain_z = NZ - (iz * (npz-1)) + 2*ghost_z;
        }
    }

	DATA_TYPE (*A)[domain_z][domain_y][domain_x] = (DATA_TYPE (*)[domain_z][domain_y][domain_x])malloc(sizeof(DATA_TYPE)*domain_x*domain_y*domain_z*2);
	if(NULL == A) return 0;

	DATA_TYPE *sbuf1 = (DATA_TYPE *)malloc(sizeof(DATA_TYPE)*domain_z*domain_y*domain_x);
	DATA_TYPE *sbuf2 = (DATA_TYPE *)malloc(sizeof(DATA_TYPE)*domain_z*domain_y*domain_x);
	DATA_TYPE *rbuf1 = (DATA_TYPE *)malloc(sizeof(DATA_TYPE)*domain_z*domain_y*domain_x);
	DATA_TYPE *rbuf2 = (DATA_TYPE *)malloc(sizeof(DATA_TYPE)*domain_z*domain_y*domain_x);

	int count1, count2;
	long traffic = 0;
	long traffic_all = 0;

#ifdef CHECK
	DATA_TYPE (*B)[domain_y][domain_x] = (DATA_TYPE (*)[domain_y][domain_x])malloc(sizeof(DATA_TYPE)*domain_x*domain_y*2);
	if(NULL == B) return 0;

	DATA_TYPE *csbuf1 = (DATA_TYPE *)malloc(sizeof(DATA_TYPE)*domain_y*domain_x);
	DATA_TYPE *csbuf2 = (DATA_TYPE *)malloc(sizeof(DATA_TYPE)*domain_y*domain_x);
	DATA_TYPE *crbuf1 = (DATA_TYPE *)malloc(sizeof(DATA_TYPE)*domain_y*domain_x);
	DATA_TYPE *crbuf2 = (DATA_TYPE *)malloc(sizeof(DATA_TYPE)*domain_y*domain_x);
#endif

	srand(100);

    init((double *)A, domain_z, domain_y, domain_x, ghost_z, domain_z - ghost_z, 
        ghost_y, domain_y - ghost_y, ghost_x, domain_x - ghost_x);

    // for (k = ghost_z; k < domain_z - ghost_z; k++) {
    //     for (j = ghost_y; j < domain_y-ghost_y; j++){
    //         for (i = ghost_x; i < domain_x-ghost_x; i++){
    //             A[0][k][j][i] = (DATA_TYPE) (1.0 * (rand() % 1024));
    //             //A[0][j][i] = 0;
    // #if  point == 0
    //             A[0][k][j][i] = ((int)A[0][k][j][i])%2;
    // #endif
    //             A[1][k][j][i] = 0;
    // #ifdef CHECK
    //             B[0][k][j][i] = A[0][k][j][i];
    //             B[1][k][j][i] = 0;
    // #endif
    //         }
    //     }
    // }


    if(tidx == 0){ //初始化x方向下边界
        for (k = ghost_z; k < domain_z - ghost_z; k++) {
            for (j = ghost_y; j < domain_y - ghost_y; j++){
                for (i = ghost_x-XSLOPE; i < ghost_x; i++){
                    A[0][k][j][i] = (DATA_TYPE) (1.0 * (rand() % 1024));
                    A[1][k][j][i] = 0;
    #ifdef CHECK
                    B[0][k][j][i] = A[0][k][j][i];
                    B[1][k][j][i] = 0;
    #endif
                }
            }
        }
	}
    if (tidx == npx - 1) { // initiate x axis bottom border
        for (k = ghost_z; k < domain_z - ghost_z; k++) {
            for (j = ghost_y; j < domain_y - ghost_y; j++){
                for (i = domain_x - ghost_x; i < domain_x - ghost_x + XSLOPE; i++){
                    A[0][k][j][i] = (DATA_TYPE) (1.0 * (rand() % 1024));
                    A[1][k][j][i] = 0;
    #ifdef CHECK
                    B[0][k][j][i] = A[0][k][j][i];
                    B[1][k][j][i] = 0;
    #endif
                }
            }
        }
    }

	if(tidy == 0){ //初始化y方向下边界
        for (k = ghost_z; k < domain_z - ghost_z; k++) {
            for (j = ghost_y-YSLOPE; j < ghost_y; j++){
                for (i = ghost_x; i < domain_x-ghost_x; i++){
                    A[0][k][j][i] = (DATA_TYPE) (1.0 * (rand() % 1024));
                    A[1][k][j][i] = 0;
    #ifdef CHECK
                    B[0][k][j][i] = A[0][k][j][i];
                    B[1][k][j][i] = 0;
    #endif
                }
            }
        }
	}
    if (tidy == npy - 1) { // initiate y axis bottom border
        for (k = ghost_z; k < domain_z - ghost_z; k++) {
            for (j = domain_y - ghost_y; j < domain_y - ghost_y + YSLOPE; j++){
                for (i = ghost_x; i < domain_x-ghost_x; i++){
                    A[0][k][j][i] = (DATA_TYPE) (1.0 * (rand() % 1024));
                    A[1][k][j][i] = 0;
    #ifdef CHECK
                    B[0][k][j][i] = A[0][k][j][i];
                    B[1][k][j][i] = 0;
    #endif
                }
            }
        }
    }

    if(tidz == 0){ //初始化z方向下边界
        for (k = ghost_z-ZSLOPE; k < ghost_z; k++) {
            for (j = ghost_y; j < domain_y - ghost_y; j++){
                for (i = ghost_x; i < domain_x-ghost_x; i++){
                    A[0][k][j][i] = (DATA_TYPE) (1.0 * (rand() % 1024));
                    A[1][k][j][i] = 0;
    #ifdef CHECK
                    B[0][k][j][i] = A[0][k][j][i];
                    B[1][k][j][i] = 0;
    #endif
                }
            }
        }
	}
    if (tidz == npz - 1) { // initiate z axis bottom border
        for (k = domain_z - ghost_z; k < domain_z - ghost_z ; k++) {
            for (j = ghost_y; j < domain_y - ghost_y; j++){
                for (i = ghost_x; i < domain_x-ghost_x; i++){
                    A[0][k][j][i] = (DATA_TYPE) (1.0 * (rand() % 1024));
                    A[1][k][j][i] = 0;
    #ifdef CHECK
                    B[0][k][j][i] = A[0][k][j][i];
                    B[1][k][j][i] = 0;
    #endif
                }
            }
        }
    }

    // 初始化各条边
    if (tidx == 0 && tidy == 0) {
        init((double *)A, domain_z, domain_y, domain_x, ghost_z, domain_z - ghost_z, 
            ghost_y - YSLOPE, ghost_y, ghost_x - XSLOPE, ghost_x);
    }
    if (tidx == 0 && tidz == 0) {
        init((double *)A, domain_z, domain_y, domain_x, ghost_z - ZSLOPE, ghost_z, 
            ghost_y, domain_y - ghost_y, ghost_x - XSLOPE, ghost_x);
    }
    if (tidy == 0 && tidz == 0) {
        init((double *)A, domain_z, domain_y, domain_x, ghost_z - ZSLOPE, ghost_z, 
            ghost_y - YSLOPE, ghost_y, ghost_x, domain_x - ghost_x);
    }
    if (tidx == npx - 1 && tidy == 0) {
        init((double *)A, domain_z, domain_y, domain_x, ghost_z, domain_z - ghost_z, 
            ghost_y - YSLOPE, ghost_y, domain_x - ghost_x, domain_x - ghost_x + XSLOPE);
    }
    if (tidx == npx - 1 && tidz == 0) {
        init((double *)A, domain_z, domain_y, domain_x, ghost_z - ZSLOPE, ghost_z, 
            ghost_y, domain_y - ghost_y, domain_x - ghost_x, domain_x - ghost_x + XSLOPE);
    }
    if (tidy == 0 && tidz == npz - 1) {
        init((double *)A, domain_z, domain_y, domain_x, domain_z - ghost_z, domain_z - ghost_z + ZSLOPE, 
            ghost_y - YSLOPE, ghost_y, ghost_x, domain_x - ghost_x);
    }
    if (tidx == npx - 1 && tidz == npz - 1) {
        init((double *)A, domain_z, domain_y, domain_x, domain_z - ghost_z, domain_z - ghost_z + ZSLOPE, 
            ghost_y, domain_y - ghost_y, domain_x - ghost_x, domain_x - ghost_x + XSLOPE);
    }
    if (tidx == 0 && tidz == npz - 1) {
        init((double *)A, domain_z, domain_y, domain_x, domain_z - ghost_z, domain_z - ghost_z + ZSLOPE, 
            ghost_y, domain_y - ghost_y, ghost_x - XSLOPE, ghost_x);
    }
    if (tidy == npy - 1 && tidz == npz - 1) {
        init((double *)A, domain_z, domain_y, domain_x, domain_z - ghost_z, domain_z - ghost_z + ZSLOPE, 
            domain_y - ghost_y, domain_y - ghost_y + YSLOPE, ghost_x, domain_x - ghost_x);
    }
    if (tidx == npx - 1 && tidy == npy - 1) {
        init((double *)A, domain_z, domain_y, domain_x, ghost_z, domain_z - ghost_z, 
            domain_y - ghost_y, domain_y - ghost_y + YSLOPE, domain_x - ghost_x, domain_x - ghost_x + XSLOPE);
    }
    if (tidy == npy - 1 && tidz == 0) {
        init((double *)A, domain_z, domain_y, domain_x, ghost_z - ZSLOPE, ghost_z, 
            domain_y - ghost_y, domain_y - ghost_y + YSLOPE, ghost_x, domain_x - ghost_x);
    }
    if (tidx == 0 && tidy == npy - 1) {
        init((double *)A, domain_z, domain_y, domain_x, ghost_z, domain_z - ghost_z, 
            domain_y - ghost_y, domain_y - ghost_y + YSLOPE, ghost_x - XSLOPE, ghost_x);
    }
    // 初始化各个角
    for (int i = 0; i < 8; i++) {
        int cnt = 0;
        int x_start, x_end, y_start, y_end, z_start, z_end;
        if (i & 0x1) { // x位的枚举符，若为0则判断tidx==0
            if (tidx == 0) {
                x_start = ghost_x - XSLOPE;
                x_end = ghost_x;
                cnt++;
            }
        } else { //非0判断tidx==npx-1
            if (tidx == npx - 1) {
                x_start = domain_x - ghost_x;
                x_end = domain_x - ghost_x + XSLOPE;
                cnt++;
            }
        }

        if ((i >> 1) & 0x1) { // y
            if (tidy == 0) {
                y_start = ghost_y - YSLOPE;
                y_end = ghost_y;
                cnt++;
            }
        } else {
            if (tidy == npy - 1) {
                y_start = domain_y - ghost_y;
                y_end = domain_y - ghost_y + YSLOPE;
                cnt++;
            }
        }

        if ((i >> 2) & 0x1) { // z
            if (tidz == 0) {
                z_start = ghost_z - ZSLOPE;
                z_end = ghost_z;
                cnt++;
            }
        } else {
            if (tidz == npz - 1) {
                z_start = domain_z - ghost_z;
                z_end = domain_z - ghost_z + ZSLOPE;
                cnt++;
            }
        }

        if (cnt == 3) { // 说明是属于某个角落的，并且start和end都设置好了
            init((double *)A, domain_z, domain_y ,domain_x, z_start, z_end, y_start, y_end, x_start, x_end);
        }
    }

	int level = 0;
	int tt,n;
	int x, y, z;
	register int ymin, ymax;
	int xmin,xmax;

    MPI_Barrier(MPI_COMM_WORLD);
    if (iam < npx * npy * npz) {
        begin = MPI_Wtime();
        for (t = 0; t < T; t++) {
            // 打包x轴方向
            count1 = pack((double*)A, domain_z, domain_y, domain_x, 
                            (tidz == 0) ? ghost_z - ZSLOPE : ghost_z, (tidz == npz - 1) ? domain_z - ghost_z + ZSLOPE : domain_z - ghost_z, 
                            (tidy == 0) ? ghost_y - YSLOPE : ghost_y, (tidy == npy - 1) ? domain_y - ghost_y + YSLOPE : domain_y - ghost_y,
                            ghost_x, 2 * ghost_x, sbuf1, t);

            count2 = pack((double*)A, domain_z, domain_y, domain_x,
                            (tidz == 0) ? ghost_z - ZSLOPE : ghost_z, (tidz == npz - 1) ? domain_z - ghost_z + ZSLOPE : domain_z - ghost_z, 
                            (tidy == 0) ? ghost_y - YSLOPE : ghost_y, (tidy == npy - 1) ? domain_y - ghost_y + YSLOPE : domain_y - ghost_y,
                            domain_x - 2 * ghost_x, domain_x - ghost_x, sbuf2, t);
            // 发送
            if (tidx & 0x1) {//奇数列进程先发送sbuf1，接收sbuf2
                MPI_Send(sbuf1, count1, MPI_DOUBLE, iam - 1, 20, MPI_COMM_WORLD);
                traffic += count1;
                MPI_Recv(rbuf2, count2, MPI_DOUBLE, iam - 1, 21, MPI_COMM_WORLD, &status);
                if (tidx != npx - 1) {
                    MPI_Recv(rbuf1, count1, MPI_DOUBLE, iam + 1, 22, MPI_COMM_WORLD, &status);
                    MPI_Send(sbuf2, count2, MPI_DOUBLE, iam + 1, 23, MPI_COMM_WORLD);
                    traffic += count2;
                }
            }
            else{//偶数列进程先接收sbuf1，发送sbuf2
                if(tidx != npx - 1){//除每行最右进程
                    MPI_Recv(rbuf1, count1, MPI_DOUBLE, iam + 1, 20, MPI_COMM_WORLD, &status);
                    MPI_Send(sbuf2, count2, MPI_DOUBLE, iam + 1, 21, MPI_COMM_WORLD);
                    traffic += count2;
                }
                if(tidx != 0){//除每行最左进程
                    MPI_Send(sbuf1, count1, MPI_DOUBLE, iam - 1, 22, MPI_COMM_WORLD);
                    traffic += count1;
                    MPI_Recv(rbuf2, count2, MPI_DOUBLE, iam - 1, 23, MPI_COMM_WORLD, &status);
                }
            }
            // 解包
            if(tidx != npx - 1){//除了每行最右进程, 解包rbuf1(横向右侧接收)
                unpack((double*)A, domain_z, domain_y, domain_x,
                        (tidz == 0) ? ghost_z - ZSLOPE : ghost_z, (tidz == npz - 1) ? domain_z - ghost_z + ZSLOPE : domain_z - ghost_z, 
                        (tidy == 0) ? ghost_y - YSLOPE : ghost_y, (tidy == npy - 1) ? domain_y - ghost_y + YSLOPE : domain_y - ghost_y,
                        domain_x - ghost_x, domain_x, rbuf1, t);
            }
            if(tidx != 0){//除了每行最左进程, 解包rbuf2(横向左侧接收)
                unpack((double*)A, domain_z, domain_y, domain_x,
                        (tidz == 0) ? ghost_z - ZSLOPE : ghost_z, (tidz == npz - 1) ? domain_z - ghost_z + ZSLOPE : domain_z - ghost_z, 
                        (tidy == 0) ? ghost_y - YSLOPE : ghost_y, (tidy == npy - 1) ? domain_y - ghost_y + YSLOPE : domain_y - ghost_y,
                        0, ghost_x, rbuf2, t);
            }

            // 打包y轴方向
            // if(tidy > 0){//除了最下层进程，打包sbuf1(纵向向下传输)
            count1 = pack((double*)A, domain_z, domain_y, domain_x, 
                        (tidz == 0) ? ghost_z - ZSLOPE : ghost_z, (tidz == npz - 1) ? domain_z - ghost_z + ZSLOPE : domain_z - ghost_z, 
                        ghost_y, 2 * ghost_y, 0, domain_x, sbuf1, t);
            // }
            // else count1 = ghost_y*domain_x;

            // if(tidy < npy - 1){//除最上层进程，打包sbuf2(纵向向上传输)
            count2 = pack((double*)A, domain_z, domain_y, domain_x, 
                        (tidz == 0) ? ghost_z - ZSLOPE : ghost_z, (tidz == npz - 1) ? domain_z - ghost_z + ZSLOPE : domain_z - ghost_z, 
                        domain_y - 2 * ghost_y, domain_y - ghost_y, 0, domain_x, sbuf2, t);
            // }
            // else count2 = ghost_y*domain_x;
            // 发送
            if(tidy & 0x1){//奇数行进程先发送csbuf1，接收csbuf2
                MPI_Send(sbuf1, count1, MPI_DOUBLE, iam-npx, 24, MPI_COMM_WORLD);
                traffic += count1;
                MPI_Recv(rbuf2, count2, MPI_DOUBLE, iam-npx, 25, MPI_COMM_WORLD, &status);
                if(tidy != npy - 1){//除每列最上进程
                    MPI_Recv(rbuf1, count1, MPI_DOUBLE, iam+npx, 26, MPI_COMM_WORLD, &status);
                    MPI_Send(sbuf2, count2, MPI_DOUBLE, iam+npx, 27, MPI_COMM_WORLD);
                    traffic += count2;
                }
            }
            else{//偶数行进程先接收csbuf1，发送csbuf2
                if(tidy != npy - 1){//除每列最上进程
                    MPI_Recv(rbuf1, count1, MPI_DOUBLE, iam+npx, 24, MPI_COMM_WORLD, &status);
                    MPI_Send(sbuf2, count2, MPI_DOUBLE, iam+npx, 25, MPI_COMM_WORLD);
                    traffic += count2;
                }
                if(tidy > 0){//除每列最下进程
                    MPI_Send(sbuf1, count1, MPI_DOUBLE, iam-npx, 26, MPI_COMM_WORLD);
                    traffic += count1;
                    MPI_Recv(rbuf2, count2, MPI_DOUBLE, iam-npx, 27, MPI_COMM_WORLD, &status);
                }
            }
            // 解包
            if(tidy != npy - 1){//除最上层进程，解包rbuf1(纵向上侧接收)
                unpack((double*)A, domain_z, domain_y, domain_x, 
                        (tidz == 0) ? ghost_z - ZSLOPE : ghost_z, (tidz == npz - 1) ? domain_z - ghost_z + ZSLOPE : domain_z - ghost_z, 
                        domain_y - ghost_y, domain_y, 0, domain_x, rbuf1, t);
            }
            if(tidy > 0){//除了最下层进程，解包rbuf2(纵向下侧接收)
                unpack((double*)A, domain_z, domain_y, domain_x, 
                        (tidz == 0) ? ghost_z - ZSLOPE : ghost_z, (tidz == npz - 1) ? domain_z - ghost_z + ZSLOPE : domain_z - ghost_z, 
                        0, ghost_y, 0, domain_x, rbuf2, t);
            }

            // z轴方向
            if(tidz > 0){//除了最下层进程，打包sbuf1(纵向向下传输)
                count1 = pack((double*)A, domain_z, domain_y, domain_x, 
                            ghost_z, 2 * ghost_z, 0, domain_y, 0, domain_x, sbuf1, t);
            }
            else count1 = ghost_z*domain_y*domain_x;

            if(tidz < npz - 1){//除最上层进程，打包sbuf2(纵向向上传输)
                count2 = pack((double*)A, domain_z, domain_y, domain_x,
                            domain_z - 2 * ghost_z, domain_z - ghost_z, 0, domain_y, 0, domain_x, sbuf2, t);
            }
            else count2 = ghost_z*domain_y*domain_x;
            // 发送
            if(tidz & 0x1){//奇数层进程先发送csbuf1，接收csbuf2
                MPI_Send(sbuf1, count1, MPI_DOUBLE, iam-npx*npy, 28, MPI_COMM_WORLD);
                traffic += count1;
                MPI_Recv(rbuf2, count2, MPI_DOUBLE, iam-npx*npy, 29, MPI_COMM_WORLD, &status);
                if(tidz != npz - 1){//除每列最上进程
                    MPI_Recv(rbuf1, count1, MPI_DOUBLE, iam+npx*npy, 30, MPI_COMM_WORLD, &status);
                    MPI_Send(sbuf2, count2, MPI_DOUBLE, iam+npx*npy, 31, MPI_COMM_WORLD);
                    traffic += count2;
                }
            }
            else{//偶数层进程先接收csbuf1，发送csbuf2
                if(tidz != npz - 1){//除最上层进程
                    MPI_Recv(rbuf1, count1, MPI_DOUBLE, iam+npx*npy, 28, MPI_COMM_WORLD, &status);
                    MPI_Send(sbuf2, count2, MPI_DOUBLE, iam+npx*npy, 29, MPI_COMM_WORLD);
                    traffic += count2;
                }
                if(tidz > 0){//除最下层进程
                    MPI_Send(sbuf1, count1, MPI_DOUBLE, iam-npx*npy, 30, MPI_COMM_WORLD);
                    traffic += count1;
                    MPI_Recv(rbuf2, count2, MPI_DOUBLE, iam-npx*npy, 31, MPI_COMM_WORLD, &status);
                }
            }
            // 解包
            if(tidz != npz - 1){//除最上层进程，解包rbuf1(纵向上侧接收)
                unpack((double*)A, domain_z, domain_y, domain_x, 
                        domain_z - ghost_z, domain_z, 0, domain_y, 0, domain_x, rbuf1, t);
            }
            if(tidz > 0){//除了最下层进程，解包rbuf2(纵向下侧接收)
                unpack((double*)A, domain_z, domain_y, domain_x, 
                        0, ghost_z, 0, domain_y, 0, domain_x, rbuf2, t);
            }

            for (z = ghost_z; z < domain_z - ghost_z; z++) {
                for (y = ghost_y; y < domain_y-ghost_y; y++){
                    for (x = ghost_x; x < domain_x-ghost_x; x++){
                        kernel(A);
                    }
                }
            }
        }
    }

    MPI_Barrier(MPI_COMM_WORLD);
    elaspe = MPI_Wtime()-begin;
	MPI_Reduce(&traffic,&traffic_all,1,MPI_LONG,MPI_SUM,0,MPI_COMM_WORLD);

    if(iam == 0)
		printf("NX = %d, NY = %d, NZ = %d, T = %d, ix = %d, ghost_x = %d, iy = %d, ghost_y = %d, iz = %d, ghost_z = %d, MStencil/s = %f, traffic = %ld\n",NX,NY,NZ,T,ix,ghost_x,iy,ghost_y,iz,ghost_z,((double)NX * NY * NZ * T) / elaspe / 1000000L,traffic_all);

#ifdef CHECK
    if(iam < npx*npy){
        for (t = 0; t < T; t++){
            count1 = 0;
            for(j = ((iam<npx)? ghost_y-YSLOPE : ghost_y); j < ((iam<npx*(npy-1)) ? domain_y-ghost_y : domain_y-ghost_y+YSLOPE); j++){//最下一行，多发送一个下边界，最上层，多发一个上边界
                for(i = ghost_x; i < ghost_x+XSLOPE; i++){
                    sbuf1[count1] = B[t%2][j][i];
                    count1++;
                }
            }

            count2 = 0;
            for(j = ((iam<npx)? ghost_y-YSLOPE : ghost_y); j < ((iam<npx*(npy-1)) ? domain_y-ghost_y : domain_y-ghost_y+YSLOPE); j++){//最下一行，多发送一个下边界，最上层，多发一个上边界
                for(i = domain_x-ghost_x-XSLOPE; i < domain_x-ghost_x; i++){
                    sbuf2[count2] = B[t%2][j][i];
                    count2++;
                }
            }

            if((iam%npx) & 0x1){//奇数列进程先发送sbuf1，接收sbuf2
                MPI_Send(sbuf1, count1, MPI_DOUBLE, iam-1, 20, MPI_COMM_WORLD);
                MPI_Recv(rbuf2, count2, MPI_DOUBLE, iam-1, 21, MPI_COMM_WORLD, &status);
                if((iam+1)%npx != 0){//除每行最右进程
                    MPI_Recv(rbuf1, count1, MPI_DOUBLE, iam+1, 22, MPI_COMM_WORLD, &status);
                    MPI_Send(sbuf2, count2, MPI_DOUBLE, iam+1, 23, MPI_COMM_WORLD);
                }
            }
            else{//偶数列进程先接收sbuf1，发送sbuf2
                if((iam+1)%npx != 0){//除每行最右进程
                    MPI_Recv(rbuf1, count1, MPI_DOUBLE, iam+1, 20, MPI_COMM_WORLD, &status);
                    MPI_Send(sbuf2, count2, MPI_DOUBLE, iam+1, 21, MPI_COMM_WORLD);
                }
                if(iam%npx != 0){//除每行最左进程
                    MPI_Send(sbuf1, count1, MPI_DOUBLE, iam-1, 22, MPI_COMM_WORLD);
                    MPI_Recv(rbuf2, count2, MPI_DOUBLE, iam-1, 23, MPI_COMM_WORLD, &status);
                }
            }
            if((iam+1)%npx != 0){//除了每行最右进程, 解包rbuf1(横向右侧接收)
                count1 = 0;
                for(j = ((iam<npx)? ghost_y-YSLOPE : ghost_y); j < ((iam<npx*(npy-1)) ? domain_y-ghost_y : domain_y-ghost_y+YSLOPE); j++){
                    for(i = domain_x-ghost_x; i < domain_x-ghost_x+XSLOPE; i++){
                        B[t%2][j][i] = rbuf1[count1];
                        count1++;
                    }
                }
            }
            if(iam%npx != 0){//除了每行最左进程, 解包rbuf2(横向左侧接收)
                count2 = 0;
                for(j = ((iam<npx)? ghost_y-YSLOPE : ghost_y); j < ((iam<npx*(npy-1)) ? domain_y-ghost_y : domain_y-ghost_y+YSLOPE); j++){
                    for(i = ghost_x-XSLOPE; i < ghost_x; i++){
                        B[t%2][j][i] = rbuf2[count2];
                        count2++;
                    }
                }
            }
            if(iam >= npx){//除了最下层进程，打包sbuf1(纵向向下传输)
                count1 = 0;
                for(j = ghost_y; j < ghost_y+YSLOPE; j++){
                    for(i = 0; i < domain_x; i++){
                        sbuf1[count1] = B[t%2][j][i];
                        count1++;
                    }
                }
            }
            else count1 = YSLOPE*domain_x;

            if(iam < npx*(npy-1)){//除最上层进程，打包sbuf2(纵向向上传输)
                count2 = 0;
                for(j = domain_y-ghost_y-YSLOPE; j < domain_y-ghost_y; j++){
                    for(i = 0; i < domain_x; i++){
                        sbuf2[count2] = B[t%2][j][i];
                        count2++;
                    }
                }
            }
            else count2 = YSLOPE*domain_x;

            if((iam/npx) & 0x1){//奇数行进程先发送csbuf1，接收csbuf2
                MPI_Send(sbuf1, count1, MPI_DOUBLE, iam-npx, 24, MPI_COMM_WORLD);
                MPI_Recv(rbuf2, count2, MPI_DOUBLE, iam-npx, 25, MPI_COMM_WORLD, &status);
                if(iam < npx*(npy-1)){//除每列最上进程
                    MPI_Recv(rbuf1, count1, MPI_DOUBLE, iam+npx, 26, MPI_COMM_WORLD, &status);
                    MPI_Send(sbuf2, count2, MPI_DOUBLE, iam+npx, 27, MPI_COMM_WORLD);
                }
            }
            else{//偶数行进程先接收csbuf1，发送csbuf2
                if(iam < npx*(npy-1)){//除每列最上进程
                    MPI_Recv(rbuf1, count1, MPI_DOUBLE, iam+npx, 24, MPI_COMM_WORLD, &status);
                    MPI_Send(sbuf2, count2, MPI_DOUBLE, iam+npx, 25, MPI_COMM_WORLD);
                }
                if(iam >= npx){//除每列最下进程
                    MPI_Send(sbuf1, count1, MPI_DOUBLE, iam-npx, 26, MPI_COMM_WORLD);
                    MPI_Recv(rbuf2, count2, MPI_DOUBLE, iam-npx, 27, MPI_COMM_WORLD, &status);
                }
            }
            if(iam < npx*(npy-1)){//除最上层进程，解包rbuf1(纵向上侧接收)
                count1 = 0;
                for(j = domain_y-ghost_y; j < domain_y-ghost_y+YSLOPE; j++){
                    for(i = 0; i < domain_x; i++){
                        B[t%2][j][i] = rbuf1[count1];
                        count1++;
                    }
                }
            }
            if(iam >= npx){//除了最下层进程，解包rbuf2(纵向下侧接收)
                count2 = 0;
                for(j = ghost_y-YSLOPE; j < ghost_y; j++){
                    for(i = 0; i < domain_x; i++){
                        B[t%2][j][i] = rbuf2[count2];
                        count2++;
                    }
                }
            }

			for (y = ghost_y; y < domain_y-ghost_y; y++){
				for (x = ghost_x; x < domain_x-ghost_x; x++){
					kernel(B);
				}
			}
        }
    }
	for (j = ghost_y; j < domain_y-ghost_y; j++){
		for (i = ghost_x; i < domain_x-ghost_x; i++){
            if(myabs(A[T%2][j][i],B[T%2][j][i]) > TOLERANCE)
				printf("%d:\tNaive[%d][%d] = %f, Check = %f: FAILED!\n", iam, j, i, B[T%2][j][i], A[T%2][j][i]);
		}
	}
#endif

    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Finalize();

}
